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Molecular dynamics simulations of 3C-SiC have been performed as a function of pressure and 
temperature. These simulations treat both electrons and atomic nuclei by quantum mechanical 
methods. While the electronic structure of the solid is described by an efficient tight-binding Hamil- 
tonian, the nuclei dynamics is treated by the path integral formulation of statistical mechanics. To 
assess the relevance of nuclear quantum effects, the results of quantum simulations are compared to 
others where either the Si nuclei, the C nuclei or both atomic nuclei are treated as classical particles. 
We find that the experimental thermal expansion of 3C-SiC is realistically reproduced by our sim- 
ulations. The calculated bulk modulus of 3C-SiC and its pressure derivative at room temperature 
show also good agreement with the available experimental data. The effect of the electron-phonon 
interaction on the direct electronic gap of 3C-SiC has been calculated as a function of temperature 
and related to results obtained for bulk diamond and Si. Comparison to available experimental data 
shows satisfactory agreement, although we observe that the employed tight-binding model tends to 
overestimate the magnitude of the electron-phonon interaction. The effect of treating the atomic 
nuclei as classical particles on the direct gap of 3C-SiC has been assessed. We find that non-linear 
quantum effects related to the atomic masses are particularly relevant at temperatures below 250 
K. 

PACS numbers: 63.20.Kr,05.30.-d,71.15.Pd 



I. INTRODUCTION 

Silicon carbide has attracted much interest in technol- 
ogy as a wide-gap semiconductor in high power electronic 
devices^ It features large electronic band gaps, extreme 
hardness, large thermal conductivity and excellent chem- 
ical stability. It is a prototype of materials which exhibit 
several polytypes (more than 200) and the only IV-IV 
compound which possesses long range order polytypes. 
The cubic modification 3C-SiC has a zinc-blende struc- 
ture, being thus the type with closest structural relation- 
ship to both diamond and elemental Si. From the point 
of view of basic research, there exists a large body of ex- 
perimental data on silicon carbide, so that the quality of 
theoretical approaches and computer simulations can be 
contrasted against them. 2 

The influence of anharmonic effects in the vibrational 
properties of 3C-SiC is revealed by several experimental 
studies. The pressure dependence of optical phonons in 
3C-SiC has been measured by means of first- and second- 
order Raman scattering up to 23 GPa, 3 while their tem- 
perature dependence has been studied by first-order Ra- 
man scattering up to 750 K.- The pressure dependence 
of phonon lifetimes, which provides detailed information 
on the anharmonic phonon-phonon interaction, has been 
also calculated from first principles using perturbation 
theory^ These calculations permitted a detailed analy- 



sis of the microscopic anharmonic mechanism and yielded 
accurate predictions of experimental properties. 

The thermal expansion of the lattice is another prop- 
erty determined by the anharmonicity of the interatomic 
potential. There is no evidence of negative thermal ex- 
pansion in 3C-SiC, in contrast to diamond and elemen- 
tal Si, which show negative thermal expansion at tem- 
peratures below 100 K2& An analysis of the most re- 
liable experimental data for the thermal expansion of 
3C-SiC led to a tabulation of recommended values for 
this quantity in Ref. [ij Theoretical results for the tem- 
perature dependence of the linear thermal expansion for 
3C-SiC have been also reported using the quasiharmonic 
approximation ^ ! 11 

The optical properties of 3C-SiC have been measured 
using spectroscopic ellipsometry as a function of tem- 
perature in the range between 90 and 550 K,^ 2 - The 
electron-phonon interaction is responsible for the de- 
crease found in the interband transition energies as tem- 
perature increases. We recall that the renormalization 
of the optical response of semiconductors by electron- 
phonon interaction has been a topic of increasing inter- 
est in recent years ) 13 i 14 with special focus in the charac- 
terization of isotopic effects^ From a theoretical point 
of view, the electron-phonon interaction in tctrahcdral 
semiconductors has been studied so far by perturbation 
theoryi&ILlS 



An interesting alternative to perturbational ap- 
proaches to study the coupling between electronic and 
vibrational degrees of freedom in solids, is the combina- 
tion of the path integral (PI) formulation with electronic 
structure methods. The path integral approach to sta- 
tistical mechanics allows us to study finite temperature 
properties that are related to the quantum nature of the 
atomic nuclei/ 19 i 20 An advantage of its combination with 
electronic structure methods is that both the electrons 
and the atomic nuclei are then treated quantum mechan- 
ically in the framework of the Born-Oppenheimcr (BO) 
approximation, so that phonon-phonon and electron- 
phonon interactions are automatically considered in the 
simulation. This unified scheme has been applied so far 
to the study of solids and molecules containing light 
atoms i 21 i 22 '2^'2£'2£'2L2£ A recent application of this 
method has shown that the electron-phonon coupling 
leads to a zero-point renormalization of the direct elec- 
tronic gap of diamond of 10%^ in agreement with a 
previous perturbational analysis^ 

In this paper we present a path integral molecular dy- 
namics study of 3C-SiC at temperatures between 100 
and 1200 K and pressures up to 60 GPa. The elec- 
tronic structure was treated with a non-orthogonal tight- 
binding (TB) Hamiltonian as a reasonable compromise to 
reduce the computational cost of deriving the BO energy 
surface for the nuclear dynamics. We are interested in the 
simulation of vibrational properties that rely on phonon- 
phonon interactions, such as the temperature dependence 
of the linear expansion coefficient, and also in the sim- 
ulation of electronic properties that are determined by 
electron-phonon coupling, such as the temperature de- 
pendence of the direct electronic gap. 

This paper is organized as follows. In Sec. II, we de- 
scribe the computational method employed in our simula- 
tions. Our results are presented and discussed in Sec. Ill, 
dealing with the thermal expansion coefficient and the 
temperature dependence of the direct electronic gap of 
3C-SiC. The results for the electronic gap in 3C-SiC will 
be then related to those derived for diamond and crys- 
talline Si. The results of the simulation will also be com- 
pared to available experimental data. The pressure de- 
pendence of the electronic gap as well as the results for 
the bulk modulus and its pressure derivative at 300 K 
will complete Sec. III. In Sec. IV, we present the main 
conclusions of the paper. 



II. COMPUTATIONAL METHOD 

The formalism employed here for the quantum treat- 
ment of electrons and nuclei is based on the combination 
of the path integral formulation, to derive properties of 
the atomic nuclei in thermal equilibrium, with an elec- 
tronic tight-binding Hamiltonian to describe the BO en- 
ergy surface, Ebo(R), of 3C-SiC as a function of the nu- 
clear configuration R. This approach has been recently 
used in our simulation of diamond 29 and isolated hydro- 



genic impurities in diamond 30 in the canonical NVT en- 
semble (number of atoms N , volume V, and temperature 
T are constant). Therefore, we present here only a brief 
summary of the method, with focus on those extensions 
required for the 3C-SiC simulations, that were performed 
in both the NVT and the isothermal-isobaric NPT en- 
semble (N, pressure P, and T are constant). The com- 
bination of the path integral formalism with ab initio 
Hamiltonians based on density funtional theory (DFT) 
is an interesting alternative that has been reviewed in 
the literature ) 31 i 32 but it has not been applied so far to 
the investigation of the temperature dependence of the 
optical response in semiconductors. Typically tight bind- 
ing methods are two orders of magnitude faster than ab 
initio DFT methods. 

The employed tight-binding one-electron effective 
Hamiltonian is based on density functional (DF) 
calculations 42. The TB energy consists of two terms, the 
first one of which is the sum of occupied one-electron 
state energies, and the second is given by a pair-wise re- 
pulsive interatomic potential. The one-electron states are 
derived by diagonalizing a two-center Hamiltonian using 
a minimal basis of non-orthogonal atomic orbitals. The 
pair potential is adjusted so that the DF energy is repro- 
duced for a series of reference systems (e.g. the dimer, the 
crystal, etc.) The total energy is thus obtained by adding 
to the electronic energy the repulsive pair-potential. Pre- 
liminary calculations on 3C-SiC revealed that the original 
parameterization of the pair-potential between Si and C 
leads to an overestimation of anharmonic effects in 3C- 
SiC, a fact that motivated us to improve the parameteri- 
zation of the Si-C pair-potential as explained in Appendix 

m 

The computational advantage of using the path- 
integral formulation of statistical mechanics is based on 
the so-called "quantum-classical" isomorphism. Thus, 
this method exploits the fact that the partition function 
of a quantum system is formally equivalent to that of a 
classical one, obtained by replacing each quantum parti- 
cle (here, atomic nucleus) by a ring polymer consisting 
of L "beads", connected by harmonic springs i 19 i 20 ' 34 i 35 
In many-body problems, the configuration space of the 
classical isomorph is usually sampled by Monte Carlo or 
molecular dynamics (MD) techniques. Here, we have em- 
ployed the PI MD method, which has been found to re- 
quire less computer time resources when applied to our 
problem. Effective algorithms to perform PI MD sim- 
ulations in the canonical NVT ensemble have been de- 
scribed in detail by Martyna et alM- and by Tuckermani^i 
The extensions for PI MD simulations in the NPT en- 
semble require the definition of appropriate dynamical 
equations of motion to include the volume as fluctuating 
dynamical variable, and multiple time step algorithms to 
integrate these equations numerically. These extensions 
are well documented in Refs. I38ll39ll40l . For the sim- 
ulation of 3C-SiC in the NPT ensemble only isotropic 
volume fluctuations were allowed. All calculations pre- 
sented here were carried out using originally developed 



software, which enables efficient PI MD simulations on 
parallel supercomputers. 

Simulations were performed on a 2 x 2 x 2 supercell of 
the 3C-SiC face-centered cubic cell with periodic bound- 
ary conditions, containing N — 64 atoms. Convergence 
criteria used previously for diamond were proved to ap- 
ply also for the 3C-SiC simulations. Thus, we used only 
the r point for the sampling of the Brillouin zone (BZ) of 
the simulation supercell in the electronic structure calcu- 
lation. A set of 4 k points would increase the computer 
time by a factor of 10 without significant changes of the 
results presented here. For a given temperature, a typical 
run consisted of 10 4 MD steps for system equilibration, 
followed by 10 5 steps for the calculation of ensemble av- 
erage properties. To have a nearly constant precision in 
the path integral results at different temperatures, we 
have taken a number of beads, L (Trotter number), that 
scales with LT = 6000 K. The atomic masses of C and Si 
were set to 12.011 and 28.086 amu, respectively. Within 
the employed formalism, the classical limit of a given 
nucleus is reached by setting the corresponding nuclear 
mass tending to infinity. Thus, for comparison with the 
results of our full PI MD simulations, we have carried out 
some PI MD simulations were the mass of either C or Si 
was set to a very large number (8 x 10 amu). The main 
effect of setting such a large nuclear mass in the path 
integral simulation is that the ring polymer associated 
to the atomic nucleus shrinks and looks just like a point 
(classical) particle. Moreover, the calculation of ther- 
mostatted equations of motion by using chains of Nose- 
Hoover thermostats (see below) ensures that the canon- 
ical probability distribution of the shrinked classical-like 
particles is correctly sampled as a function of tempera- 
ture. Also classical MD simulations were performed with 
the same interatomic interaction. The classical limit is 
easily achieved within the PI algorithm by setting the 
Trotter number L = 1. 

The quantum simulations were performed using a stag- 
ing transformation for the bead coordinates. Chains of 
four Nose-Hoover thermostats were coupled to each of 
the staging variables to generate the canonical NVT 
ensemble^ To integrate the equations of motion we have 
used the reversible reference system propagator algo- 
rithm (RESPA), which allows one to define different time 
steps for the integration of the fast and slow degrees of 
freedom^ For the evolution of the fast dynamical vari- 
ables, that include the thermostats and harmonic bead 
interactions, we used a time step St = At/A, where At is 
the time step associated to the calculation of TB forces. 
A value of At = 0.5 fs was found to provide adequate 
convergence. The thermostat "mass" parameter, Q, is 
chosen to evolve in the scale of the harmonic bead forces 
by being defined as^ 
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FIG. 1: Temperature dependence of the cell parameter of 
3C-SiC. Open circles are the results of PI MD simulations at 
zero pressure. The continuous line is a fit to Eq. ((3}- The 
dotted line is the result derived from the experimental values 
of the thermal expansion coefficients given in Ref. |9|. The 
error bar is derived from the uncertainty of the experimental 
data. The statistical error of the simulation results is of the 
size of the symbols. 



was coupled to the volume and the barostat "mass" pa- 
rameter was set 100 times larger than the thermostat 
one. By calculating the average pressure in the NVT en- 
semble we could check the internal consistency between 
our NPT and NVT simulations. The pressure estimator 
used in the NVT ensemble is given in Appendix [Bj The 
direct electronic gap of 3C-SiC was derived as 



Ea — (E c ) — (E v 



(2) 



where (3 = (ksT) 1 represents the inverse temperature. 
In the case of the NPT ensemble a chain of four barostats 



where (E c ) and (E v ) are the ensemble average of the one- 
electron states associated with the bottom of the con- 
duction band and the top of the valence band at the T 
(k = 0) reciprocal lattice point. For details on the calcu- 
lation of the expectation values (E c ) and (E v ) see Ref. 

Several sets of simulations were performed in this 
work. Two sets of 3C-SiC simulations were done at 
temperatures between 100 and 1200 K. One set corre- 
sponds to NPT simulations at constant pressure (P=0), 
and the other to NVT simulations at constant volume 
(a = 4.3594 A). Another set of NPT simulations was 
made at room temperature (T ~ 300 K) and pressures in 
the range between -10 and 60 GPa. Additional constant 
volume simulations were done for diamond and Si as a 
function of temperature. 



III. RESULTS AND DISCUSSION 

The results derived from our PI MD simulations are 
presented in the next subsections. Firstly, we focus on 
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FIG. 2: The temperature dependence of the cell parameter of 
3C-SiC (circles) derived from PI MD simulations is compared 
to PI MD simulations performed by selectively setting either 
the Si nuclei (squares) or the C nuclei (diamonds) as classical 
particles. The results of classical MD simulations are shown 
by triangles. The lines represent empirical fits as explained 
in the text. The statistical error of the simulation results is 
of the size of the symbols. 

the temperature dependence of the cell parameter and 
the linear expansion coefficient at zero pressure. 



A. Thermal expansion 

The temperature dependence of the cell parameter, 
a(T), of 3C-SiC at zero pressure is presented in Fig. Q] 
The simulation results are plotted as open circles. The 
dotted line is derived from the recommended experimen- 
tal values of the linear expansion coefficient of 3C-SiC by 
a numerical integration.— The resulting relative changes 
in the cell parameter, Aa(T), were converted to an abso- 
lute scale by using the experimental value of a — 4.3596 
A at 297 K.— Both simulation and experimental data 
were fitted with a standard Bosc-Einstein expression^ 



a{T) = a + b 1 + 



exp(9 a /T) - 1 



(3) 



the values of the fitted parameters are summarized in 
Table[U The results in Fig. [T]show that the simulations 
and experimental data agree within the uncertainty of the 
experimental values. At temperatures below 600 K the 
calculated cell parameter a(T) shows a rigid shift with 
respect to experiment of about 10 -3 A toward larger val- 
ues. This shift increases slightly at higher temperatures. 
The most plausible explanation for this high temperature 
deviation of simulated data from experiment is that the 
employed TB model overestimates the anharmonicity of 
the interatomic potential when the amplitude of nuclei 



TABLE I: Values of the parameters obtained by fitting the 
temperature dependence of the lattice constant, a(T), of 3C- 
SiC with Eq. ((3j . The results correspond to NPT simulations 
at P = 0. The experimental parameters were derived from 
the recommended thermal expansion coefficients of Ref. @. 
The PI MD parameters for those simulations where the Si 
or C nuclei are treated in the classical limit were obtained 
by adding a linear term, aiT, to Eq. @. The classical MD 
parameters correspond to the cubic function given in Eq. Q . 





do (A) 


oi (Ak- 1 ) 


6(A) 


e«(K) 


PI MD 

exp. 

PI MD (Si cla.) 

PI MD (C cla.) 


4.3463 
4.3475 
4.3458 
4.3451 


1.0 10" 5 

1.1 10" 5 


0.0131 
0.0110 
0.0097 
0.0077 


975.7 
898.0 
1077.2 
921.7 




a (A) 


oi (Ak- 1 ) 


a 2 (Ak- 2 ) 


a 3 (Ak- 3 ) 


classical MD 


4.3483 


2.1 10~ 5 


6.4 10~ 9 


-2.3 1Q- 12 



vibrations is activated thermally to values significantly 
larger than low temperature zero-point motions. 

It is interesting to quantify the renormalization of the 
cell parameter a(T) as a consequence of the quantum 
character of the atomic nuclei by comparison to results 
of classical simulations. Moreover, the employed formal- 
ism allows us to treat the atomic nuclei either quantum 
mechanically or in the classical limit just by tuning the 
nuclear mass. Thus, we have also calculated the effect 
in a(T) of considering only one type of nuclei classically. 
The results of these simulations are presented in Fig. [5] 
The full PI MD results (circles) are compared to classi- 
cal MD simulations (triangles). While the classical cal- 
culations yield a finite linear expansion for T tending 
to zero, the quantum PI MD calculations yield cell pa- 
rameters which become independent of T at low T, in 
agreement with the experimental results. The quantum 
PI MD results that selectively treat the Si or the C nu- 
clei as classical particles are shown by squares and dia- 
monds, respectively. The lines represent numerical fits to 
the simulation results. The classical MD data were fitted 
with a cubic function, 



(l-cl 



«cn = 



E 



a,T 



(4) 



while the PI MD simulation results treating either the Si 
or C nuclei classically were fitted to Eq. $5§ plus an addi- 
tional linear term, a\T, to account for the finite slope of 
the curves at T = 0. The values of the fitted coefficients 
are summarized in Table [I] 

The cell parameter a c i a (T = 0) amounts to 4.3483 
A while the extrapolated value of the PI MD simula- 
tion is a(T = 0)= 4.3594 A. Therefore the zero-point 
renormalization of the cell parameter of 3C-SiC amounts 
to Aa/a — 2.5 x 10~ 3 (Aa = a — a c / a )- This value is 
slightly lower than the average (2.9 x 10 -3 ) of the zero- 
point renormalization of diamond (3.9 x 10~ 3 ) and Si 
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FIG. 3: Temperature dependence of the cell parameter renor- 
malization of 3C-SiC. The results derived from full quantum 
simulation (continuous line) are compared to those derived 
by treating either the Si or the C nuclei in the classical limit 
(dashed and dotted lines, respectively). The dashed-dotted 
line shows the sum of the separate Si and C renormalizations. 
The curves have been derived from the fits shown in Fig. [2] 



(1.9 x 10 3 ) derived from experiments on crystals with 
different isotopic masses J£ 

The a(T) curves where either the C or Si nuclei are 
treated classically provide evidence for the linear char- 
acter of the renormalization of the lattice parameter. 
This fact is clearly seen in Fig. [3j where the absolute 
values of the cell parameter renormalization, Aa(T) = 
a(T) — a c i a (T), are represented up to 800 K. We observe 
that the sum of the cell parameter renormalizations ob- 
tained when either the Si or C nuclei are treated quan- 
tum mechanically is nearly identical to the total renor- 
malization obtained in the full PI MD simulations. This 
linear behavior is probably related to the fact that the 
cell parameter renormalization is relatively small, so that 
in terms of perturbation theory it can be realistically 
described by second-order terms in the atomic displace- 
ments. 

The recommended experimental values of the linear 
thermal expansion, a(T), of 3C-SiC are shown in Fig. 
2] as closed circles. Error bars represent the uncertainty 
of the experimental results. The simulation results pre- 
sented in Fig. [4] were derived from the fits shown in 
Fig. [2] for a(T). The agreement between our PI MD re- 
sults for a(T) and the experimental data is satisfactory 
in the whole temperature range. The largest deviation 
is found at temperatures above 700 K where the thermal 
expansion of 3C-SiC is overestimated by about 8% by 
our model. Although the deviation between PI MD sim- 
ulated results and measured data is well within the exper- 
imental error bar, we stress that the TB model seems to 
overestimate the anharmonicity of the interatomic poten- 
tial at temperatures above 600 K. The thermal expansion 






& 
X 



a 2- 



- 




1 ' 


.._ 


_ T 


- 


SiC, P = 


















--* • • J-^*"*^ > 1 






+ • 4- jT • 


__ _ ■ — ■ 






/ / • „•'-" J 




/ /• . • " 




■ 1 ¥ K 


. / / _£. 


■ I/I 


: I It 


••' ' i 


: I 4 1 


I J ■ 


I I 








_ • 


classical MD 


_ 




PIMD(C classical) 






PI MD (Si classical) 




j 






" 




J 


• exp. 


- 


II 


DFT 




til 






yj i , i 



"0 500 1000 

temperature (K) 

FIG. 4: Thermal expansion coefficient of 3C-SiC as a function 
of temperature. Closed circles are experimental data from 
Ref. [SJ and their estimated error bars are shown at 500 and 
1000 K. The simulation results are derived from the numerical 
fits shown in Fig. [2] The classical MD result at temperatures 
above 800 K requires to add higher-order terms to the cubic fit 
given by Eq. Q. The results derived from DFT calculations 
in combination with a quasiharmonic approximation are taken 
from Ref. [U 



coefficient obtained in the classical MD simulations at a 
given temperature is always larger than the correspond- 
ing PI MD value. This behavior is expected from the 
fact that the classical thermal expansion is finite in the 
zero-temperature limit (while it vanishes in the quantum 
case) and from the consideration that both sets of simu- 
lations should converge one to the other at high enough 
temperatures. Note that the deviation from experiment 
is always larger for classical than for quantum simula- 
tions, while the PI MD results where either the Si or the 
C nuclei are treated as classical particles lie between both 
limits. The classical simulation results presented in Fig. 
[4] have required to add higher-order terms to the cubic 
fit given by Eq. (j4|) at temperatures above 800 K. 

The quality of our results is comparable to previous 
calculations of a(T) by using either ab initio DFT elec- 
tronic Hamiltonians 1 * or phenomenological models^ to 
determine harmonic vibrational frequencies and a quasi- 
harmonic approximation to take into account the anhar- 
monicity of the interatomic potential. The DFT result of 
Ref. 11 has been plotted in Fig. 4. The deviation from 
experimental data found at high temperature is probably 
caused by anharmonic effects not taken into account by 
the quasiharmonic approximation. 



B. Direct electronic gap 

The one-electron energies derived by the TB model for 
3C-SiC at the symmetry point T are compared to results 



TABLE II: Relative one-electron energies at the symmetry 
point r for a static 3C-SiC lattice. LMTO results arc from 
Ref. |H, and EPM data from Ref. 0. All values in eV. 
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of linear muffin-tin-orbital calculations^ 2 - and empirical 
pseudopotential method^ in Table Hfl These energies are 
obtained with the atoms fixed in their crystallographic 
positions and thus neglect the effects of lattice vibra- 
tions. The employed TB model predicts that the first 
direct gap, E Q , of 3C-SiC at T appears between elec- 
tronic states with Tis symmetry. This fact is in contra- 
diction with the other electronic Hamiltonians that show 
that the conduction band bottom is of Ti symmetry at 
r, 42 i 43 In the case of diamond and Si the bottom of the 
conduction band at T is found to be of Tis symmetry 
by ab initio calculations! 44 ' 45 The repulsion between the 
T25' valence band and the Tis conduction band, induced 
by the asymmetric potential, is probably responsible for 
the lowering of the Ti conduction band with respect to 
the Ti5 counterpart in SiC. Although the employed TB 
model does not provide an accurate description of the 
conduction bands, it has demonstrated, by the study of 
the direct gap of diamond, that is a realistic starting 
model to study electron-phonon interaction effects^ 

The temperature dependence of the direct gap, E (T), 
derived by our NPT simulations of 3C-SiC is shown by 
circles in Fig. [5] The continuous line represents a fit of 
Eq (T) to the Bose-Einstein expression 



E (T) = e -9 1 + 



exp(e B /T) - 1 



(5) 



The extrapolated value of Eq(T — 0) amounts to 6.9 eV, 
while at 300 K the gap is reduced to 6.84 eV, and at 
1000 K amounts to 6 eV. The decrease of Eq with tem- 
perature is an effect of the electron-phonon interaction, 
and therefore it is expected to depend on the amplitude 
of the nuclei displacements. To assess this point, the re- 
sults obtained for E^{T) by treating cither the Si or the 
C nuclei as classical particles are presented in Fig. [5] as 
squares and diamonds, respectively. In both cases the 
results were fitted to Eq. $5§ with an additional linear 
term, eiT, to account for the finite slope at T = 0. The 
fitted parameters are summarized in Table IIIII The re- 
sults of classical MD simulation for Eo tC i a (T) are given 
by triangles in Fig. O while the dashed-dotted line is a 
cubic fit as in Eq. (fj|. The classical simulation predicts 
Eo,da{T = 0) = 7.59 eV, thus the zero-point renormal- 
ization of the direct energy gap, E — E ^ c [ a , amounts to 
-0.69 eV, which is roughly 10% of the value of the gap. 
The zero-point renormalization in 3C-SiC is similar to 
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FIG. 5: Temperature dependence of the direct energy gap, 
E (T), of 3C-SiC at T at P = 0. The results derived from full 
quantum simulation (circles) are compared to PI MD simula- 
tions performed by selectively setting the Si nuclei (squares) 
or the C nuclei (diamonds) as classical particles. The results 
of classical MD simulations are shown by triangles. The lines 
represent empirical fits as explained in the text. The statisti- 
cal error of the simulation results is less than the size of the 
symbols. 



TABLE III: Values of the parameters obtained by fitting 
the temperature dependence of the direct energy gap, Eq(T), 
of 3C-SiC with Eq. ©. All results were derived by NPT 
simulations at P = 0. The PI MD parameters for those cases 
where the Si or C nuclei are treated in the classical limit 
were obtained by adding a linear term, e\T, to Eq. |[5j|. The 
classical MD parameters correspond to a cubic function as in 
Eq. ©. 





eo (eV) 


ei (eV K" 1 ) 


ff(eV) 


e*(K) 


PIMD 

PI MD (Si cla.) 
PI MD (C cla.) 


7.548 
7.578 
7.646 


-2.7 10~ 4 
-3.6 10" 4 


0.643 
0.501 
0.363 


887.2 
817.9 
615.7 




e (eV) 


ei (eV K" 1 ) 


e 2 (eV K" 2 ) 


e 3 (eV K" 3 ) 


classical MD 


7.593 


-12.1 10~ 4 


-6.2 10~ 7 


3.2 10~ 10 



that found for diamond by either PI MD simulations 29 
or by perturbation theory^ 

The temperature dependence of the direct gap renor- 
malization, AE a (T) = E Q (T) - E , c i a (T) is shown in Fig. 
H] by a full line. The separate Si contribution to the 
gap renormalization is obtained by a simulation where 
the C nuclei are treated classically. This result is shown 
by a dotted line in Fig. [6l The dashed line represents 
the separate C contribution as derived from a simulation 
with classical Si nuclei. The sum of both Si and C incre- 
ments (dashed-dotted line) is somewhat larger than the 
gap renormalization obtained in the full quantum simu- 
lation of 3C-SiC, in particular at temperatures below 250 
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FIG. 6: Temperature dependence of the direct band gap 
renormalization, AEo(T), of 3C-SiC at P = 0. The results 
derived from full quantum simulation (continuous line) are 
compared to those derived by treating either the Si or the C 
nuclei in the classical limit (dashed and dotted lines, respec- 
tively). The dashed-dotted line shows the sum of the separate 
Si and C renormalizations. The curves have been derived from 
the fits shown in Fig. [5] 



FIG. 7: Temperature dependence of the direct energy gap, 
Eo(T), of 3C-SiC at V at constant volume. The results derived 
from full quantum simulation (circles) are compared to PI 
MD simulations performed by selectively setting either the 
Si nuclei (squares) or the C nuclei (diamonds) as classical 
particles. The results of classical MD simulations are shown 
by triangles. The lines represent empirical fits as explained 
in the text. The statistical error of the simulation results is 
less than the size of the symbols. 



TABLE IV: Values of the parameters obtained by fitting 
the temperature dependence of the direct energy gap, Eq(T), 
of 3C-SiC. The results correspond to constant volume NVT 
simulations using a cell parameter of a — 4.3594 A. See Tab. 
IIIII for details on the fitting functions. 





e (eV) 


ei (eV K" 1 ) 


fl(eV) 


6 E (K) 


PI MD 

PI MD (Si cla.) 

PI MD (C cla.) 


7.461 
7.442 
7.508 


-2.3 10 -4 
-3.2 10" 4 


0.555 
0.379 
0.252 


876.4 
722.6 
517.5 




eo (eV) 


ei (eV K- 1 ) 


e 2 (eV K" 2 ) 


e 3 (eV K" 3 ) 


classical MD 


7.541 


-11.8 10~ 4 


-4.2 10~ 7 


2.4 10~ 10 



K, where the gap renormalization reaches its largest val- 
ues. This non-linear behavior of AEq (T) is in contrast to 
the linear one found for the lattice parameter renormal- 
ization, Aa(T), in Fig. [3l This non-linearity is probably 
related to the fact that the relative gap renormalizations 
are found to be much larger than the relative cell param- 
eter renormalizations. In terms of perturbation theory, 
the fact that the total gap renormalization is lower that 
the separate Si and C contributions, implies that fourth- 
order terms in the nuclei displacements are important for 
AEq(T), and their contribution is of opposite sign to the 
leading second-order terms. 

The conclusion that fourth-order terms in the gap 
renormalization are not related to changes in the cell pa- 
rameter is further demonstrated by the results shown in 
Figs. [7]and[51 where we show the values of the direct gap, 
Eq(T), and gap renormalizations, AEq(T), derived at 



constant volume by NVT simulations. The volume was 
kept fixed at the equilibrium value of the cell parameter 
(a = 4.3594 A) extrapolated for T = from our PI MD 
simulations (see Fig. [1]). The direct gap, Eo(T), in Fig. 
[7] is presented for the full PI MD simulation, the classical 
MD simulation, and PI MD simulations treating either 
the Si or the C nuclei as classical particles. For each case 
numerical fits to the simulation results were performed in 
the same way as explained for the zero-pressure results in 
Fig. [5] The fitted parameters are summarized in Table 
IIVI The gap renormalizations, AEo(T), in Fig. [8] show 
that non-linear effects appear again at temperatures be- 
low 250 K in the case of constant volume simulations. 
Note that the comparison between Figs. [6] and [8] reveals 
that non-linear effects in the gap renormalization are in- 
dependent of the volume. 

The difference between the values of Eo (T) obtained in 
our PI MD simulations at P = (Fig. [5]) and at constant 
volume (Fig. [7]) allows us to quantify the thermal expan- 
sion effect in the gap. This difference is shown in Fig. [9] 
as a function of temperature. We note that the thermal 
expansion produces a decrease in the value of Eo(T) as 
temperature increases. At 300 K this decrease amounts 
to only 8 meV, while at 1000 K, the decrease rises to a 
value slightly less than 110 meV. These values should be 
compared to the total effect of temperature in the direct 
gap, E a (T) - E a (T = 0), which amounts to 64 meV at 
300 K and to 900 meV at 1000 K, as derived from the PI 
MD results shown in Fig. \5\ Thus, the thermal expan- 
sion is responsible for about 10 % of the total decrease 
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FIG. 8: Temperature dependence of the direct band gap 
renormalization, AEo(T), of 3C-SiC at constant volume. The 
results derived from full quantum simulation (continuous line) 
are compared to those derived by treating either the Si or 
the C nuclei in the classical limit (dashed and dotted lines, 
respectively). The dashed-dotted line shows the sum of the 
separate Si and C renormalizations. The curves have been 
derived from the fits shown in Fig. [7] 
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FIG. 9: Thermal expansion effect in the value of the direct 
electronic gap, Eo, of 3C-SiC as a function of temperature. 
The circles are the difference between the Eo(T) values ob- 
tained by PI MD simulations at P — (see Fig. [5} and at 
constant volume (see Fig. 0. The line is a guide to the eye. 



in the value of E at temperatures of 300 and 1000 K. 



C. Comparison of Eo for diamond, 3C-SiC and Si 

It is interesting to compare the simulation results of 
the direct gap Eq(T) found for 3C-SiC with those cor- 
responding to diamond and Si by using the same tight- 
binding parameterization. The studied direct gap, Eq, 
for both diamond and Si corresponds to transitions be- 
tween one-electron states with symmetry 1^25' (valence 
band) and F^ (conduction band). To prevent possible 
inaccuracies of the TB model in the determination of the 
thermal expansion of diamond and Si, we have performed 
constant volume simulations of diamond and Si with the 
following values of the cell parameters, ac = 3.567 A for 
diamond^ and asi = 5.430 A for Si4£ Our simulation 
results of Eq(T) in diamond, 3C-SiC, and Si will be also 
compared to available experimental data. Note that in 
this comparison the thermal expansion effect in Eq(T) is 
not included in our constant volume simulations. 

In Fig. [10] we show the results of the PI MD and clas- 
sical MD simulations for the relative shifts of the direct 
gap of diamond with temperature. The lines through 
the points correspond to numerical fits using Eqs. (|5"|) 
and [[4"]). respectively. The fitted parameters are summa- 
rized in Table [Vj The zero-point renormalization of Eq 
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FIG. 10: Relative shifts of the direct electronic gap of dia- 
mond obtained by our PI MD (circles) and classical MD sim- 
ulations (triangles) as a function of temperature. The dashed 
line is a fit to the experimental data of Logothetidis el al. for 
diamond IIa4£ The zero of the energy scale was set at 7.06 
eV for both experimental and simulation results. 



TABLE V: Values of the parameters obtained by fitting the 
temperature dependence of the direct energy gap, Eo(T), of 
diamond, 3C-SiC, and Si with Eq. ©. The PI MD results 
correspond to constant volume NVT simulations. 





eo (eV) 


ff(eV) 


e B (K) 


PI MD C (a=3.567 A) 


7.802 


0.740 


1665.0 


PI MD 3C-SiC (a=4.3594 A) 


7.461 


0.555 


876.4 


PI MD Si (a=5.4296 A) 


3.511 


0.204 


513.6 


exp. C a 


7.387 


0.320 


1060.0 


exp. 3C-SiC 6 


7.943 


0.230 


668.0 


exp. Si c 


3.467 


0.091 


382.6 


exp. Si d 


3.378 


0.025 


267.0 



"Logothetidis et al, Ref. |4g| 
6 Petalas et al, Ref. Q3 
c Jcllison and Modine, Ref. |5CI 
d Lautenschlager et al, Ref. 49 



amounts to 0.705 eV, a value that agrees well with the 
result of 0.678 eV derived by Zollner et al. in Ref. |T7| by 
a perturbational treatment of the electron-phonon cou- 
pling. Unfortunately there is no experimental estimation 
of the renormalization of the direct gap, Eq, of diamond. 
However, the zero-point renormalization of the indirect 
gap of diamond derived from luminescence data amounts 
to 0.37 eV.— The broken line shown in Fig. [10] is the 
experimental result reported in Ref. 48 for diamond Ha, 
based on measurements of the complex dielectric function 
by spectroscopic ellipsometry between 100 to 650 K. The 
extrapolated experimental value Eq(T = 0) varies from 
7.06 to 7.14 eV, depending on the line-shape analysis of 
the spectra by using a first or a second derivative. The 
TB model gives a value of E (T — 0) — 7.06 eV in good 
agreement with experiment. The slope of the PI MD 
results at temperatures above 500 K is larger than that 
of the experimental data, a fact that indicates that the 
employed TB model overestimates the electron-phonon 
interaction at temperatures above 500 K. 

In Fig. QT]we compare the relative shifts of the direct 
energy gap of 3C-SiC obtained by our PI MD and clas- 
sical MD simulations (circles and triangles, respectively) 
to a fit to the experimental data derived between 90 and 
550 K by spectroscopic ellipsometry. The experimental 
data can not discriminate between intcrband electronic 
transitions occurring at the points T (critical point Eq), 
and along A (critical point E\) in reciprocal space, as 
they appear in the same energy region. The fitted pa- 
rameters are summarized in Table |V] The comparison to 
experiment shows that our computational model tends 
to overestimate the shift in the energy gap. 

The simulation results for the temperature shift of the 
direct electronic gap of Si are compared to available ex- 
perimental data in Fig. [T2j The dashed line represents 
the numerical fit to the spectroscopic ellipsometric data 
of Lautenschlager et a/.j 4 ^ while the dotted line is de- 
rived by polarization modulation ellipsometry!^ The dif- 
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FIG. 11: Relative shifts of the direct electronic gap of 3C-SiC 
obtained by our PI MD (circles) and classical MD simulations 
(triangles) as a function of temperature. The dashed line is a 
fit to the experimental data of Petalas et al^ The zero of the 
energy scale corresponds to 6.91 eV for the simulation results 
and to 7.61 eV for the experimental ones. 



ference between both sets of experimental data might be 
due to the fact that the measured excitations are a su- 
perposition of intcrband transitions (Eo, E\) along the 
A direction that includes both the T and L points at its 
boundary. Fitted parameters are collected in Table [V] 
The extrapolated PI MD value of E (T = 0) amounts 
to 3.31 eV, in reasonable agreement to the experimen- 
tal extrapolated results of 3.35 eV, 49 and 3.38 eV^ The 
calculated zero-point renormalization of Eq amounts to 
0.15 eV. This result is to be compared to the experimen- 
tal value of 0.12±0.02 eV derived by Lastras-Martmez et 
al. from a study of isotopically pure and natural Si^ 
This experimental value has to be considered as a mix- 
ture of Eq and E\ transitions. The temperature shifts in 
Eq derived from our PI MD simulations appear again to 
be larger than in the experimental data, pointing toward 
an overestimation of the electron-phonon interaction by 
the employed electronic TB Hamiltonian. 

In order to compare the effects of the electron-phonon 
interaction in the direct gap, Eq, of diamond, 3C-SiC, 
and Si, we have plotted in Fig. [T3]the simulation results 
of Eq(T) using relative energy and temperature scales. 
For each crystal the energy was measured in units of its 
direct gap Eq(T = 0) and the temperature in units of the 
corresponding Debye temperature. We see that by using 
reduced units the direct gap of 3C-SiC falls roughly be- 
tween the values found for diamond and Si in the stud- 
ied temperature range. It has been recently shown that 
at very low temperatures (below 10 K for Si) the gap 
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FIG. 12: Relative shifts of the direct electronic gap of Si 
obtained by our PI MD (circles) and classical MD simulations 
(triangles) as a function of temperature. The dashed line is a 
fit to the experimental data of Lautenschlager ei a/.,— while 
the dotted line is a fit to the Jellison and Modine results 
based on polarization modulation ellipsometry.— The zero of 
the energy scale is set at 3.31 eV for the simulation results, 
at 3.35 eV for the dashed line and at 3.38 eV for the dotted 
line. 
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changes with T like T 4 . 52 Our calculations, however, do 
not have the necessary accuracy to reveal this depen- 
dence. 



D. Pressure dependence of Eq in 3C-SiC 

The pressure dependence of the direct gap of 3C-SiC at 
300 K was derived from simulations in the NPT ensem- 
ble up to 60 GPa. The results of the quantum PI MD and 
classical simulations are shown in Fig. Q3] The main dif- 
ference between both sets of results is a rigid shift of the 
Eq values, that reflects the dependence of the electron- 
phonon coupling on the nuclei displacements around the 
equilibrium positions. At a given pressure, the vibra- 
tional amplitudes are always larger for the quantum sim- 
ulations. The cell parameter difference, a versus a c i a , has 
a relatively smaller effect in the shift of Eq. For example, 
at P = the quantum result for the direct gap is 0.35 eV 
lower than the classical one (see Fig. ITU) , a value that 
includes the effect of the volume difference in the quan- 
tum and classical simulations. The corresponding result 
obtained at a constant volume at 300 K is 0.31 eV (see 
Fig. 0. 

The derivative, dEg/dP, at P = is readily obtained 
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FIG. 14: Pressure dependence of the direct electronic gap of 
3C-SiC as derived from PI MD and classical MD simulations 
in the NPT ensemble at 300 K. The statistical error of the 
simulation results is less than the size of the symbols. 
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FIG. 15: P-V curves derived for 3C-SiC PI MD and classical 
MD simulations in the NPT ensemble at 300 K. The lines are 
fits to the Murnaghan's equation of state using the points in 
the interval between [-10,10] GPa. The statistical error of the 
simulation results is of the size of the symbols. 



from the data in Fig. [21 giving a result of 39 meV/GPa 
in the quantum case versus a value of 41 meV/GPa in the 
classical limit at 300 K. These values are in reasonable 
agreement with the calculation of Park et aL,— which 
gives a value of 51 meV/GPa, based on ab initio cal- 
culations using a local-density-functional approximation 
(LDA) without considering any kind of temperature ef- 
fect. 

Finally, the set of simulations performed at 300 K al- 
lows us to plot P — V curves that can be used to de- 
rive the value of the bulk modulus, -Bo, and its pressure 
derivative, B' , for 3C-SiC. In Fig. OH we show the P-V 
curves obtained from NPT simulations of 3C-SiC for sev- 
eral pressures in the range between -10 and 30 GPa. The 
five points calculated in the pressure interval [-10,10] GPa 
were used to fit a Murnaghan equation of stated 



-DO 



-1/Bo 



(6) 



where the subindex indicates values at P = 0. The 
continuous line in Fig. [15] shows the fit to the PI MC 
simulations, that provides the values B =225 GPa and 
B' = 4.1. We checked that this value of Bq is consistent 
with that (Bq = 221 ± 8 GPa) derived by calculating the 
volume fluctuations in the NPT ensemble at P = 0. The 
fluctuation relation for a given pressure is 



B = k B T 



(V) 



(v 2 ) - {vy 



(7) 



The calculated value of Bo is in excellent agreement 
with the experimental data ) 55 i 56 ' 57 and also the calcu- 
lated value of B' Q = 4.1 agrees with experiment (_B = 



4 ± 0.3) j££ The values obtained from the Murnaghan fit 
of the classical simulation results were B c ; a =230 GPa 
and B' = 4.2. 



IV. CONCLUSIONS 

The simulation method employed in this work has 
demonstrated its capability for the description of anhar- 
monic effects related to the phonon-phonon interaction, 
as well as for the treatment of the electron-phonon cou- 
pling in semiconducting solids as a function of temper- 
ature. Thus, this type of simulations is an alternative 
to perturbational treatments, with the advantage of be- 
ing also applicable in cases where the convergence of a 
perturbational series might be slow. A prerequisite to 
account for phonon-phonon and electron-phonon inter- 
actions is a quantum description that includes both elec- 
trons and nuclei. In this respect, the Feynman formu- 
lation of statistical mechanics allows us to simulate the 
quantum mechanical properties of the atomic nuclei in 
the solid at finite temperatures with either Monte Carlo 
or Molecular Dynamics techniques. A particular advan- 
tage of PI MD is that the algorithm can be easily paral- 
lelized, as the calculation of total energies and forces for 
each of the L different replicas of the solid, that are gen- 
erated by the discretization of the path integral (Trotter 
number), can be run in a different processor. The use 
of an electronic Hamiltonian to describe the interatomic 
interactions in the solid allows us to study also electronic 
properties in the simulations. At present we have limited 
ourselves to simplified, although accurate, Hamiltonians 
of the tight-binding type. But in the future, an interest- 
ing improvement will be the combination of PI with ab 
initio parameter-free Hamiltonians. 

The PI MD simulation of 3C-SiC has been able to 
reproduce different experimental properties of the solid. 
Moreover, the comparison to classical simulations allowed 
us to evaluate the magnitude of quantum effects as a 
function of temperature, as well as to determine the sep- 
arate contribution of each type of nuclei (either C or 
Si) to several properties. The employed potential model 
predicts a cell parameter as a function of temperature, 
a(T), in good agreement with the available experimental 
data. The correct description of this anharmonic prop- 
erty points toward a realistic treatment of the phonon- 
phonon interactions by our computational model. The 
main discrepancy is found at temperatures above 600 K, 
where the linear expansion coefficient is predicted to be 
about 8% larger than the experimental one. The zero- 
point renormalization of the cell parameter of 3C-SiC is 
calculated to be Aa/a = 2.5 x 10~ 3 (Aa = a — a c i a ). 
We find that this value is close to the sum of zero-point 
renormalizations obtained when only one type of atomic 
nuclei (either Si or C) is treated quantum mechanically. 
This near linear behavior is probably related to the fact 
that the cell parameter renormalization depends only on 
second-order terms in the nuclei displacements, due to 
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the small magnitude of Aa/a. 

The direct electronic gap at T of 3C-SiC has a signif- 
icant temperature dependence as a consequence of the 
electron-phonon coupling. The effect of the zero-point 
vibrations of the lattice phonons leads to a gap renor- 
malization of A^o = -0.69 eV (AE = E - E , c i a ). 
This effect is so large that any theoretical approach aim- 
ing at a quantitative determination of the direct elec- 
tronic gap can not be based only on an improved solu- 
tion of the many-body electron problem, but it should 
also include the treatment of the electron-phonon inter- 
action. The calculated relative value of the zero-point 
renormalization is AEq/Eq = 0.10. In this case, the 
sum of the separate Si and C contributions is found to 
be lower than the total zero-point renormalization. This 
non-linear behavior suggests that fourth-order terms in 
the nuclei displacements are important in this case, as 
a result of the large value of AEq/Eq. The ratio of the 
partial contributions of C and Si to the zero-point gap 
renormalization is 1.67 at P = 0. This is rather close to 
the inverse square root of the ratio of the corresponding 
masses (28/12) 1 / 2 = 1.53. The latter represents the ra- 
tio of squares of zero-point vibrational amplitudes, under 
the assumption of equal force constants, in the harmonic 
approximation. The calculated pressure coefficient of Eq 
is 39 meV/GPa at 300 K. An ab initio calculation with- 
out considering temperature and electron-phonon inter- 
actions gives a value of 51 meV/GPa. Our calculation 
of the bulk modulus of 3C-SiC at 300 K (B = 225 GPa) 
and its pressure derivative (B' = 4.1) shows quantitative 
agreement with experiment. 

The calculated zero-point renormalization of the direct 
gap at r of diamond and Si amounts to 0.7 eV and 0.15 
eV, respectively. For Si the renormalization derived from 
spectroscopic ellipsometry of isotopic crystals amounts to 
0.12±0.02 eV. 51 The experiment can not discriminate be- 
tween direct transitions at T and along A, as they appear 
at similar energies. In the case of diamond, a calculation 
based on perturbation theory results in a value of the 
direct gap renormalization of 0.68 eV^ in good agree- 
ment with our non-perturbational result. The experi- 
mental value for the renormalization of the indirect gap 
of diamond amounts to 0.37 eV^ The comparison of the 
simulation results of diamond, 3C-SiC and Si with avail- 
able experimental data shows that the employed model 
Hamiltonian tends to overestimate the decrease of Eq 
with temperature, i.e., the electron-phonon coupling in 
the employed TB model seems to be too strong. This 
fact is another motivation for improving the electronic 
model in future work. 
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APPENDIX A: PAIR POTENTIAL FOR SI-C 

The pair-potential is parameterized by using the fol- 
lowing functional form 33 



HI 



E SlC (R) = }_^d n (R c -Ry 

n=2 



(Al) 



where R is the interatomic Si-C distance, R C =AA au, and 
Esic(R) = if R > R c . The original coefficients d n were 
modified to decrease the anharmonicity of the Si-C po- 
tential at distances around 3.55 au, that corresponds to 
the nearest- neighbors in 3C-SiC. The modified d n coef- 
ficients are (in au): d 2 = 0.06825, d 3 = -0.49329, d 4 = 
2.37716, d 5 = -5.79511, d 6 = 7.90779, d 7 = -6.27467, 
d 8 = 2.86625, d 9 = -0.69579, d 1Q = 0.06926. The pro- 
cedure to fix the d n coefficients was to calculate curves 
of internal energy versus volume in the classical limit at 
T = 0, i.e., the atoms occupy fixed equilibrium positions 
at a given volume. The Taylor expansion of the internal 
energy around its minimum at V r j a can be expressed as 

A£~I Bda /Am2 



2 Vda 



(AVf + --^(-l-B> cla )(AV)\ (A2) 



6 VI 



with AV = V — V c ia, and B c i a , and B' cla being the bulk 
modulus and its pressure derivative in the classical limit 
at T = and P — 0. The values used in this expansion 
were B cla = 245 GPa, B' cla = 4, and a c i a = V^ = 4.3457 
A. The original analytic form of Esic(R) was modified 
to obtain a reasonable approximation to the AE curve 
in Eq. (|A2[) for a range of volumes defined by a cell 
parameter in the interval a c i a ± 0.4 au. 



APPENDIX B: PRESSURE ESTIMATOR IN THE 
NVT ENSEMBLE 

The cartesian coordinates of the N atomic nuclei in the 
crystal are denoted as x\ a , where the superscript a runs 
from 1 to 3 A. The subindex i denotes the "bead" asso- 
ciated to a given atomic nucleus and runs from 1 to the 
Trotter number L. The staging coordinates u" are de- 
fined by a linear transformation of x\ a that diagonalizes 



the harmonic energy between neighboring bead; 



,(«) _ Joe) 



,37,39 



(Bl) 



(a) (a) {i ~ 1) (a) 1 (a) 

a i — x i ■ x i+l ■ x l 



, * 



,L. (B2) 
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The pressure estimator in our NVT PI MD simulation 
is obtained as 



1 



L 3N 



"-&XE- 



» 



W 



,(«) 



3V m r,2^ 



1 yy dU(R t 



» 



,(«) 



(B3) 

where mj™' 1 and v 4 w represent the dynamic mass and 

velocity associated to the staging coordinate u\ a . These 
masses are given a o 37 i 39 

m^ = m a , (B4) 



i- 1 



^f^Ct •> ^ ^5 . . , • ±J 



(B5) 



where m a is the nuclear mass associated to coordinate 
a. E arm represents the harmonic energy between beads, 



that in terms of the stagging coordinates is written as 



L 3N 



^arm. — 



JED 

i=2 a=l 



» 



,(«) 



(B6) 



The last summand in Eq. (|B3[) represents the volume 
derivative of the potential energy of the crystal for the 
nuclei configuration R; = (xj, . . . ,x? N ). We have per- 
formed this derivative numerically by considering an ex- 
pansion of 10 -4 A in the cell parameter, 2a, of the sim- 
ulation cell, although it is also possible to calculate this 
derivative from the stress tensor, obtained through the 
Hellman-Feynman theorem. The fractional coordinates, 
xf 1 /2a, of the nuclei remain constant along this volume 
expansion, implying that all cartesian coordinates, xf , 
must change according to the modified value of 2a. 
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